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PREFACE 


The  investigation  reported  herein  was  conducted  by  personnel  of  the 
Geomechanics  Division  (GD),  Structures  Laboratory  (SL) ,  U.  S.  Army  Engineer 
Waterways  Experiment  Station  (WES) .  It  was  sponsored  by  the  Defense  Nuclear 
Agency  under  Task  Y99QAXSB,  "Ground  Shock  Predictions,"  Work  Unit  00020, 
"Waveform  Comparison  Techniques." 

The  study  was  conducted  and  this  report  prepared  and  written  by 
Dr.  G.  Y.  Baladi  and  Mr.  D.  E.  Barnes  (GD)  during  the  period  October  1981- 
October  1982  under  the  general  direction  of  Mr.  Bryant  Mather,  Chief,  SL,  and 
Dr.  J.  G.  Jackson,  Jr.,  Chief,  GD. 

COL  Tilford  C.  Creel,  CE,  was  Commander  and  Director  of  WES  during  the 
investigation  and  publication  of  this  report.  Mr.  F.  R.  Brown  was  Technical 
Director. 
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CONVERSION  FACTORS,  METRIC  (SI)  TO  U.  S.  CUSTOMARY 
UNITS  OF  MEASUREMENT 
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U.  S.  customary  units  as  follows: 
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CHAPTER  1 


INTRODUCTION 


1.1  BACKGROUND 

It  has  been  and  still  is  customary  to  analyze  explosion-generated  ground 
shock  waveforms  (measured,  computed,  or  both)  subjectively.  This  is  accom¬ 
plished  by  comparing  two  or  more  waveforms  and  verbalizing  their  compati¬ 
bility  through  the  aid  of  statements  such  as  "the  peaks  are  within  a  factor 
of  two"  or  "the  overall  agreement  is  pretty  good."  Each  analyst,  however, 
has  his  own  opinion  about  what  "pretty  good"  may  mean,  and  their  opinions 
quite  often  differ  greatly.  Consequently,  subjective  discrepancy  measures 
have  probably  produced  as  much  confusion  and  controversy  as  they  have 
enlightenment  on  a  host  of  ground  shock  issues. 

It  is  time  to  minimize  the  confusion  and  controversy.  Waveforms  should 
be  compared  objectively,  using  discrepancy  measures  that  are  rooted  in 
statistical  theory.  This  report  treats  two  such  measures  and  recommends  them 
for  adoption  by  the  ground  shock  calculation/measurement  community. 

Within  the  framework  of  the  theory  of  probability  there  are  two 
approaches  that  can  be  taken  to  develop  objective  waveform  discrepancy 
measures.  The  first  approach  involves  straightforward  application  of  statis¬ 
tical  concepts  to  obtain  ensemble  average  (mean),  mean  square,  standard 
deviation,  etc.,  for  a  given  instant  of  time.  To  use  this  approach,  however, 
it  is  necessary  to  have  information  about  the  probability  distribution  of  a 
ground  shock  parameter  throughout  the  time  history  of  its  response  or  at  least 
a  large  number  of  individual  responses  or  measurements  obtained  at  the  same 
location.  The  second  approach  involves  the  use  of  temporal  averages  and 


temporal  mean  squares  in  order  to  compare  two  response  histories  and  make  an 
objective  judgement  on  their  agreement  or  disagreement  throughout  a  given 
period  of  time  or  "time  window." 

Using  the  second  approach,  T.  L.  Geers  (Reference  1)  developed  two 
objective  discrepancy  measures  for  comparing  transient  response  histories; 
these  were  the  temporal  root  mean  square  and  the  correlation  error  history 
measure.  The  objective  discrepancy  measures  developed  in  this  report  closely 
parallel  Geers'  development. 

1.2  OBJECTIVE 

The  primary  objective  of  this  study  was  to  develop  and  document  obje 
tive  waveform  discrepancy  measures  for  comparing  arbitrary  transient  res; 
histories.  Secondary  objectives  were  (a)  to  incorporate  the  newly-developed 
waveform  discrepancy  measures  into  a  computer  program  which  can  read  digi¬ 
tized  measured  or  calculated  waveforms  and  produce  objective  waveform  compari 
sons  and  perform  probabilistic  analyses  on  a  given  number  of  response  time 
histories,  and  (b)  to  demonstrate  the  potential  utility  of  the  computer 
program  using  the  results  of  recent  field  experiments  and  code  calculations. 

1.3  SCOPE 

The  theoretical  development  behind  statistical  objective  discrepancy 
measures  is  presented  in  Chapter  2.  Chapter  3  demonstrates  the  application 
of  the  objective  discrepancy  measures  through  the  use  of  simple  analytic 
sinusoidal  waveforms.  To  demonstrate  the  capabilities  of  the  computer 
program  WCT  (Waveforms  Comparison  Technique),  statistical  analyses  of 
measured  data  and  examples  of  how  calculated  response  histories  can  be 
compared  to  measurements  are  given  in  Chapter  4.  Chapter  5  summarizes  the 
report  and  presents  recommendations. 


CHAPTER  2 


STATISTICAL  METHOD  FOR  COMPARISON  OF  TRANSIENT 
RESPONSE  HISTORIES 


2.1  INTRODUCTION 

In  general,  a  waveform  is  characterized  by  its  amplitude  and  its  fre¬ 
quency.  Thus,  the  comparison  of  two  waveforms  must  be  approached  with  these 
features  in  mind.  In  addition,  phase  shifts  must  be  considered. 

Historically,  the  shock  and  vibrations  community  has  characterized 
individual  waveforms  by  assigning  them  an  average  amplitude  and  by  decom¬ 
posing  their  frequency  content  to  obtain  a  mean  square  spectral  density 
function  (References  2,  3,  and  4).  The  average  amplitude  most  commonly 
employed  has  been  the  root  mean  square  value.  Similar  concepts  and  parame¬ 
ters  are  used  in  the  following  sections  to  develop  objective  waveform  dis¬ 
crepancy  measures. 

2.2  BASIC  EQUATIONS 
2.2.1  Single  Waveform 

Let  P(t)  be  a  periodic  function  of  period  T  .  Under  very  general 
conditions,  P(t)  may  be  represented  by  a  superposition  of  sinusoids  using 
the  following  exponential  Fourier  series  (Reference  5) : 


P(t)  =  SieXP  ^nWo^  (2.J 

n=-°° 

where  i  =  /-I  is  a  complex  number,  wq  =  2tt/T  is  the  fundamental  angular 
frequency,  and  Cn  is  Fourier  coefficients  that  can  be  evaluated  directly 
from  the  relation 

1  T//-2 

C  =  —  /  P(t)exp  (-inw  t)  dt  (2.. 


Using  Equation  2.1  and  Parseval's  theorem  (Reference  5),  it  can  be  shown  that 


1 

T 


P2(t)  dt  - 


E  N2 

ns-co 


(2.3) 


Note  that  the  left-hand  side  of  Equation  2.3,  called  the  temporal  mean  square 
of  P(t)  ,  equals  the  sum  of  the  squares  of  the  absolute  values  of  the 
Fourier  coefficients.  Hence,  the  temporal  mean  square  is  indicative  of  the 
amplitude  of  P(t)  . 

2.2.2  Two  Waveforms 

Let  P^(t)  and  +  ^  two  identical  waveforms  except  for  a  con¬ 

stant  phase  shift  between  them  (equal  to  <f>  ) .  Such  waveforms  are  correlated 
Reference  3  defines  this  correlation  as  the  temporal  autocorrelation  function 
X(<f>),  where 


T/2 

X(*>  -  Y  J  Pl(t)P2(t  +  dt  (2.4) 

-T/2 

Note  that  when  4>  ■  0  ,  P^(t)  *  P£(t)  =  P(t)  ,  and  Equation  2.4  reduces  to 

the  temporal  mean  square  of  P(t)  . 

Because  x  C«J> )  is  related  to  the  mean  square  spectral  density  function 
(Reference  3)  which  determines  the  frequency  decomposition  of  a  given  wave¬ 
form,  Equation  2.4  is  indicative  of  the  frequency  content  of  waveforms  as 
well  as  their  phase  shifts. 

2.3  OBJECTIVE  DISCREPANCY  MEASURES 

Based  on  Equations  2.3  and  2.4,  T.  L.  Geers  (Reference  6)  suggested 
three  objective  discrepancy  measures  for  comparing  two  (numerically  or 
experimentally  generated)  waveforms. 


Consider  R^(t)  to  be  an  errorless  or  true  response  function  and  I^Ct) 
to  be  a  similar  response  history,  but  they  differ  somewhat  in  amplitude, 
frequency,  and  phasing.  Geers  defined  two  correlation  factors  to  charac¬ 
terize  the  differences  between  and  R2  in  terms  of  (a)  magnitude  (i.e., 

amplitude),  and  (b)  phase  and  frequency;  namely. 


(2.5) 


and 


Pcf(t) 


f t  R1(t)R 
Jo 

2(t)  dr 

(t)  dr 

1/2 

X  r2(t)  dT. 

1/2 

(2.6) 


Here,  Mc^(t)  is  the  magnitude  correlation  factor,  and  P^f(t)  is  the 
phase -and -frequency  correlation  factor.  Note  the  distinct  preservation  of 
the  above  fundamental  character  of  Equations  2.3  and  2.4  in  the  above 
expressions. 

Geers  also  defined  a  combined  correlation  factor  to  enfold  the  magni¬ 
tude  correlation  factor  and  the  phase-and-f requency  correlation  factor  into 
one  expression,  i.e.. 


cef(t) 


(t)  - 


(2.7: 


Finally,  the  magnitude  error,  phase-and-frequency  error,  and  combined  error 
were  defined  by  Geers  as 


(2.8) 
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mag 


(t) 


E 


phs 


(t) 


E  (t)  =  SIGN 
com 


-  Mcf(t>  -  1 
'  1  -  Pcf(t) 

{[W'>]2  +  [v(t)]2}1/2 


(2.9) 

(2.10) 


Equations  2.8  through  2.10  represent ’powerful  measures  for  quantifying 
temporal  discrepancies  between  given  waveforms;  however,  because  they  all 
involve  time  integrations,  they  are  discrepancy  measures  throughout  a  given 
time  window  rather  than  time-discrete  measures.  This  offers  certain  advan¬ 
tages  because  the  quality  of  waveforms  throughout  their  time  histories  is 
what  is  important  in  designing  a  structure  to  sustain  such  waveforms. 

Note  that  the  definition  of  E  (t)  in  Equation  2.10  capitalizes  on 
the  orthogonality  of  (t)  and  ®-pjls(t^  »  as  shown  in  Figure  2.1  and 

defined  by  Equations  2.8  and  2.9.  Also,  in  keeping  with  Figure  2.1,  it  can 
be  easily  shown  (using  Equations  2.5  and  2.6)  that 


0  < 


EPhs(t) 


(2.11) 


2.4  ENSEMBLE  AVERAGING 

Quite  often,  situations  arise  in  which  several  waveforms  need  to  be 
compared  as  a  group;  e.g.,  when  redundant  field  records  and/or  multiple 
calculations  are  available.  The  above  error  concepts  can  readily  be  extended 
to  cover  these  situations  by  "ensemble  averaging." 

For  N  records  in  a  set  (either  calculated  or  measured),  average  or 
mean  error  factors  may  be  defined  as 


12 


^  •->  *-•-  '■L 


^ m 


MEAN 


(2.: 


A  great  advantage  occurs  in  using  Equation  2.14  (rather  than  straight¬ 
forward  statistical  methods)  to  compute  the  mean  combined  error;  i.e.,  one 
avoids  the  calculation  of  standard  deviations  (and  other  statistical  mea¬ 
sures)  for  E  (t)  and  E  ,  (t)  .  This  is  due  to  the  vector  magnitude 
mag  phs 

aspect  of  E  (t)  ;  i.e.,  ±E  (t)  or  ±E  ,  (t)  produces  the  same 

com  mag  phs 


E 

com 


(t)  . 


In  Chapter  3  we  demonstrate  the  utility  of  Equations  2.8  through  2.10 
and  Equations  2.12  through  2.14  by  applying  them  to  analyses  of  simple 


sinusoidal  waveforms. 
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CHAPTER  3 


ANALYTIC  EXPOSITION  OF  OBJECTIVE  DISCREPANCY  MEASURES 


3.1  INTRODUCTION 

In  this  chapter  the  statistical  measures  described  in  Chapter  2  (Equa¬ 
tions  2.8  through  2.10)  are  examined  analytically  using  three  pairs  of  con¬ 
trived  sinusoidal  waveforms.  In  Section  3.2,  two  undamped  waveforms  are  used 
to  demonstrate  the  objective  description  of  phase  and  magnitude  discrepan¬ 
cies;  Section  3.3  extends  this  analysis  to  include  a  frequency  discrepancy. 
Section  3.4  adds  the  further  complication  of  slight  damping. 

3.2  EXAMPLE  1;  UNDAMPED  SINUSOIDAL  RESPONSE; 

PHASE  AND  MAGNITUDE  DIFFERENCES 

Consider  the  following  two  sinusoidal  responses  (Figure  3.1): 


R^(t)  *  sin  2irt 


R-(t)  -  (1  +  e  )  sin  (2irt  +  <|>) 
z  m 


(3.1) 


(3.2) 


where  t  is  time  in  milliseconds.  Assume  that  R^(t)  is  an  errorless  base 
or  true  response  while  R2(t)  is  a  comparable  response  history  with  an  error 
in  magnitude  equal  to  ,  and  an  error  in  phase  equal  to  <f>  .  Substitution 
of  Equations  3.1  and  3.2  into  Equations  2.8  and  2.9  leads  to 


sin  2irt 


E  (t) 
mag 


sin  2irt 


cos  (2irt  +  2«f>) 

I  Tl72~ 


1  +  e 


(3.3) 


COS  2TTt 


sin  2irt 


sin  2wt 


cos  (2iTt  +  2<fc) 


cos  (2irt  +  <{>) 


sin  2irt 


cos  2?rt 


SINUSOIDAL  RESPONSE 


The  combined  error  can  be  calculated  directly  using  Equation  2.10  and  the 
results  of  Equations  3.3  and  3.4. 


Note  that  for  large  values  of  t  (t  >  2  in  this  problem).  Equations 
3.3  and  3.4  rapidly  approach  limits,  i.e.,  they  become 


and 


1  +  e 


m 


1 


(3.5) 


COS  d> 


(3.6) 


respectively.  Consequently,  Equation  2.10  also  approaches  a  limit.  These 

limits  are  indicated  on  Figures  3.2  through  3.5  which  illustrate  the  behavior 

of  Equations  2.10,  3.3,  and  3.4  for  this  example  (in  which  e  =  0.5  and 

m 

$  *  0.6  radian).  It  is  clear  from  these  figures  that  within  a  very  short 
time  the  objective  discrepancy  measures  have  essentially  captured  the  correct 
values  of  the  magnitude  and  phase  errors  and  therefore  improve  their  acqui¬ 
sition  with  time. 

As  a  final  note,  if  <J>  =  0  and  e  >_  -1  ,  Equations  3.3  and  3.4  (as 
well  as  Equations  3.5  and  3.6)  reduce  to 


and 


e 
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(3.7) 
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Moreover,  for  t  >  2  ,  Equation  3.5  can  be  rewritten  as 


E 

mag 


for 


E  > 

m  — 


-1 


1 


E _ (t)  =  -(e_  +  2) 


for  e  <  -1 


(3.9) 


which  leads  to 


_1  i  Emag(t)  i  "  (3-10> 

Further,  for  t  >  2  ,  Equations  3.6  and  3.8  indicate  that 

0-EPh.<Oil  <3-U) 

which  is  a  conclusion  that  was  previously  stated  in  Equation  2.11. 

And,  finally,  note  that  if  the  absolute  value  brackets  were  to  be 
omitted  from  the  numerator  of  the  fraction  in  Equation  2.6,  the  present 
example  problem  would  yield 


1  +  e 

E  ,  (t)  1 nr— — “T"  cos  <j> 

phs  1  +  e 

m 


(3.12) 


which  would  make  the  phase  error  dependent  upon  .  This,  in  turn,  could 

lead  to  unreliable  results.  For  example,  if  d>  =  0  and  e  =  -1  +  6  , 

m 

where  6  is  a  small  positive  increment  <<1  ,  Equation  3.12  gives 

E  ,  (t)  ~  0  ;  yet  for  <p  =  0  and  e  =  -1  -  6  ,  E  .  (t)  ~  2  .  This 
phs  m  phs 

•  Ephs(t)  cal- 
culations  (without  the  absolute  value)  might  be  unreliable. 


suggests  that  in  practical  cases  with 


R2(t) 


<< 


R]_(t) 


3.3  EXAMPLE  2;  UNDAMPED  SINUSOIDAL  RESPONSE; 

ERASE,  FREQUENCY  AND  MAGNITUDE  DIFFERENCES 

In  this  example,  the  following  sinusoidal  responses  are  considered 


(Figure  3.5) 


I 


R^(t)  =  sin  2tt t 


(3.13) 
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=  (1  +  em)  sin  |^2rr ( 1  +  e^)  t  +  4>J 


(3.14) 


Like  the  previous  example,  R^(t)  is  assumed  to  be  an  errorless  base  or  true 

response  while  R^Ct)  is  a  comparable  response  history  with  error  in  magnitude 

equal  to  e  ,  error  in  frequency  equal  to  e,  ,  and  error  in  phase  equal  to 
‘  m  i 

<J>  .  For  this  example,  e  =  0.5,  <p  =  0.6,  and  er  =  0.005. 

m  t 

Substitution  of  Equations  3.13  and  3.14  into  Equations  2.8  and  2.9 
leads  to 


E(t)  = 
mag 


sinj27t(l  +  e^)t 
1  «  /  ^  .  \ 


n  j  att  ^  x  ■+*  z.)  t  ip 

—  C0SL27T(1  +  £f)I:  +  2*. 
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in  27tt 
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1/2 


1  +  e  -  1  (3.15) 
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2irt 


cos  2irt 


and 


Ephs(t)  ■  1  - 
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sin  2neft 

sin  2tt (2  + 

-  sin  $ 

r  2 

sin  n  e^t 

2  - 
sin  ii  (2  +  cf)t 

2ttef  t 

2  it  (2  +  ef)t 

t 

n(2  +  6f)t 

(3.16) 


,  sin  2irt 

1  ’  -ST"  cos 


2tf  t  J 


i/2 


2n(l 


and,  as  before,  the  combined  error  can  be  calculated  using  Equation  2.10 

and  the  results  of  Equations  3.15  and  3.16.  Figures  3.6  through  3.8 

present  the  behavior  of  Equations  2.8  through  2.10  for  this  example  (in 

which  e  =0.5,  <$>  =  0.6,  and  tc  =  0.005). 

m  t 

For  t  >  2  and  <<  1  ,  Equations  3.15  and  3.16  reduce  to  (see 


Figures  3.6  through  3.8) 


TIME,  MSEC 


(3.18) 


sin  7t £. £ t 

- - —  cos  (4>  +  ire..t) 

7re,t  T  f 


Note  that  Equation  3.17  is  identical  to  the  corresponding  result  for  the 
previous  example  (Equation  3.5).  Note  too  that  if  2  <  t  < < ( tt e ^ )  ^  , 

Equation  3.18  is  essentially  identical  to  its  earlier  counterpart 
(Equation  3.6);  i.e.,  Ep^s(t)  ~  1  -  jcos  <fr|;  however,  for  t  >>  (e^)  ^ 

Ephs  ~  1  (see  Figure  3.7). 

From  this  example  we  conclude  that  frequency  error,  as  embodied  in 

the  term  ,  has  a  negligible,  if  any,  effect  on  E  (t)  ,  yet  has  a 

i  mag 

profound  effect  on  E^^t)  •  In  order  to  put  this  into  proper  context, 
however,  the  effects  of  damping  must  be  considered. 

3.4  EXAMPLE  3;  LIGHTLY  DAMPED  SINUSOIDAL  RESPONSE; 

PHASE,  FREQUENCY  AND  MAGNITUDE  DIFFERENCES 

In  this  example  the  effects  of  damping  on  the  sinusoidal  responses  of 
the  second  example  (i.e..  Equations  3.13  and  3.14)  are  considered.  We 
rewrite  Equations  3.13  and  3.14  (with  damping)  as  (Figure  3.9) 

R^(t)  =  sin  (2irt)exp  (-8t)  (3.19) 

and 

R2(t)  *  (1  +  tm)  sin  [2tt(1  ■+•  e^)t  +  <|>]exp  [-(1  +  £d)8t]  (3.20) 

where  e  ,  4>  ,  and  er  are  as  before,  8  is  the  damping  factor  (0.4 

m  r 

msec  1  in  this  example)  and  is  the  error  in  damping  (assumed  to  be  0.1 

for  this  example).  The  effects  of  the  damping  in  Equations  3.19  and  3.20  can 
be  clearly  seen  by  comparing  Figures  3.9  and  3.5. 

Substitution  of  Equations  3.19  and  3.20  into  Equations  2.8  and  2.9  leads 
to 
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LEGEND 


- R^(t)  =  sin(2irt)  exp(-gt) 


Figure  3.9 


Example  3:  Time  histories  of  two  lightly  damped 
sinusoidal  responses;  %  =  0.5  ,  <J>  =  0.6  radian 


0.005 


0.1,  and  0  *  0.4  (msec)~l. 


The  combined  error  (Equation  2.10)  becomes 


E 

com 


(t) 


/k-  /a 
j /k  -  f&\ 


/ab 


(3.27) 


Figures  3.10  through  3.12  present  the  behaviors  of  Equations  3.21,  3.22,  and 
3.27. 

For  t-*»  ,  Equations  3.21  and  3.22  become 


Equations  3.21  and  3.22  or  Equations  3.28  and  3.29  indicate  that  both  the 
magnitude  error  and  the  phase-and-f requency  error  are  dependent  on  the  damping 
parameter  8/2*  .  However,  the  phase-and-frequency  error  is  independent  of 


sot.  that  If  c2  «  1  ,  e2  «  1  ,  1  ,  *2  «  1  ,  6  <  1  ,  and 


(2*ef/B)  <  1  ,  Equations  3.28  and  3.29  become 
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Figure  3.11  Time  history  of  phase- and- frequency  error  for 
example  3. 
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8  sin  2<j>  •L/z 

L  2tt  J 

(3.31) 


In  this  case,  the  magnitude  error  depends  only  on  e  and  e,  ;  and  the 

m  a 

phase-and-frequency  error  depends  on  ,  <fi  ,  and  the  damping  parameter 
0/2ir  . 

The  previous  three  examples  demonstrated  the  application  of  the  objec¬ 
tive  discrepancy  measures.  The  potential  utility  of  the  computer  program  WCT 
is  demonstrated  in  the  next  chapter  through  statistical  analyses  of  measured 
data  and  examples  of  how  calculated  response  histories  can  be  compared  to 


measurements . 


CHAPTER  4 


APPLICATIONS 


4.1  INTRODUCTION 

The  objective  discrepancy  measures  established  in  Chapter  2  (Equations 
2.5  through  2.13)  were  incorporated  into  a  computer  program  named  WCT.  The 
listing  of  WCT,  its  flow  chart,  and  its  user's  guide  are  included  in  Appendix 
A.  The  computer  program  WCT  is  capable  of  processing  digitized  data  tapes 
containing  either  measured  or  calculated  waveforms  to  produce  (a)  the  mean 
value  and  standard  deviation  at  each  time-step  of  any  set  of  transient 
response  histories,  and  (b)  time  histories  of  the  objective  discrepancy 
measures  established  in  Chapter  2  (Equations  2.5  through  2.13)  for  any  pair 
of  waveforms.  To  demonstrate  this  capability,  WCT  was  used  to  analyze 
selected  free-field  data  recorded  on  the  DISC  Test  I  and  II  events  (References 
7  and  8),  which  were  High  Explosive  Simulation  Technique  (HEST)  experiments 
performed  in  the  desert  alluvium  of  Ralston  Valley,  Nevada. 

4.2  STATISTICAL  ANALYSIS  OF  MEASURED  DATA 

Figure  4.1  presents  nine  cavity  pressure  measurements  and  their  inte¬ 
grals  recorded  for  DIS  Test  I.  These  cavity  pressure  measurements  were 
input  to  the  WCT  code  which  integrated  them  and  produced  the  mean  integral 
and  its  standard  deviation  bounds,  as  shown  in  Figure  4.2.  The  mean  integral 
was  then  differentiated  to  obtain  the  mean  cavity  pressure  waveform  (also 
shown  in  Figure  4.2).  This  mean  cavity  pressure-time  history  and  its 
standard  deviation  bounds  were  subsequently  used  as  airblast  pressure  drivers 
for  one-dimensional  ground  shock  calculations  of  the  DISC  Test  II  event,  as 
reported  in  Reference  9. 
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Early-time  histories  of  mean  cavity  pressure  and  integral  with  standard 
deviation  bounds  for  integrals;  DISC  Test  I  event. 


4.3  COMPARISON  OF  MEASURED  VERSUS 

CALCULATED  RESPONSE  HISTORIES 

Using  the  statistical  variations  of  cavity  pressure  and  soil  compressi¬ 
bility  from  DISC  Test  I  as  input,  a  series  of  probabilistic  ID  ground  shock 
calculations  was  performed  to  predict  particle  velocity  at  the  3-meter  depth 
for  the  DISC  Test  II  event  (Reference  9) .  The  expected  value  obtained  from 
these  calculations  and  three  records  of  measured  velocities  are  plotted  in 
Figure  4.3.  Subjectively  the  comparison  looks  "pretty  good."  But  in  order 
to  obtain  a  more  objective  judgment  of  the  degree  of  agreement  or  disagree¬ 
ment  among  these  velocities,  the  waveforms  of  Figure  4.3  were  input  to  the 
computer  program  WCT,  using  the  expected  value  from  the  probabilistic  ID 
calculations  as  a  base  (truth)  record.  The  time  histories  of  the  magnitude 
errors,  the  phase-and-frequency  errors,  and  the  combined  errors  are  shown  in 
Figures  4.4  through  4.6,  respectively.  The  errors  associated  with  two  of  the 
measured  waveforms  (the  dotted  and  the  dash-dotted  curves)  are  uniformly 
small  and  essentially  identical.  The  errors  are  larger  for  the  dashed  curve, 
but  only  during  the  initial  2-  to  3-msec  toe  (or  precursor) .  The  ensemble 
averages  of  the  individual  errors  in  Figures  4.4  through  4.6  are  shown  in 
Figures  4.7  through  4.9.  These  are  simply  the  mean  values  of  the  errors 
computed  using  Equations  2.12,  2.13,  and  2.14.  Comparison  of  Figures  4.7  and 
4.8  indicates  that  the  dominant  errors  in  this  case  are  the  magnitude  errors. 

Figures  4.10  through  4.13  compare  the  mean  of  the  three  DISC  Test  II 
measurements  (dashed  curve)  with  the  calculated  expected  value  (solid  line) . 
The  magnitude  error  has  a  plus-and-minus  oscillation  during  the  rise  portion 
and  then  settles  on  a  numerical  value  of  minus  0.1.  For  all  practical 
purposes,  the  phase-and-frequency  error  is  essentially  zero;  consequently, 
the  combined  error  is  dominated  by  the  magnitude  error. 
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Figure  4.3  Comparison  of  calculated  and  measured  pa*'  le  velocity-time 
histories  for  DISC  Test  II  event;  depth  =  3.0  metres. 
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Figure  4.4  Time  histories  of  magnitude  errors  between  measured 
and  calculated  velocity  waveforms;  DISC  Test  II  event, 
depth  *  3.0  metres. 
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Figure  4.5  Time  histories  of  phase-and-f requency  errors  between 
measured  and  calculated  velocity  waveforms;  DISC 
Test  II  event,  depth  *  3.0  metres. 
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-  GAGE  NO.  2307,  RANGE  =6.0  METRES 
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Figure  4.6  Time  histories  of  combined  error  (magnitude,  phase  and 
frequency)  for  each  measured  velocity  waveform  relative 
to  the  calculated  waveform;  DISC  Test  II  event,  depth  = 
3.0  metres. 
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Figure  4.7  Time  history  of  magnitude  error  between  mean  measured 
and  calculated  velocity  waveforms;  DISC  Test  II  event 
depth  =  3.0  metres. 
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Fisure  4.8  Time  history  o£  phase-and-f requency  error  b«““n 
Figure  ‘>asured  calcuiated  velocity  waveforms;  BISC  Test 

event,  depth  *=  3.0  metres. 


TIME,  MSEC 


p,olir.  A  Q  Time  history  of  combined  error  (magnitude  and  phase-and- 
Figure  4.9  Time^histo  y^^  ^  neasured  and  calculated  velocity 

waveforms;  DISC  Test  II  event,  depth  -  3.0  metres. 
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CHAPTER  5 


SUMMARY  AND  RECOMMENDATIONS 

5.1  SUMMARY 

A  set  of  objective  discrepancy  measures  for  the  comparison  of  transient 
response  histories  has  been  established.  It  consists  of  the  magnitude 
correlation  factor,  the  phase-and-frequency  correlation  factor,  the  magnitude 
error,  the  phase-and-frequency  error,  and  the  combined  magnitude  and  phase-and- 
frequency  errors.  Their  validity  and  behavior  were  checked  and  demonstrated 
for  several  simple  sinusoidal  responses. 

The  objective  discrepancy  measures  were  incorporated  into  a  computer 
program  named  WCT  which  processes  digitized  data  tapes  containing  measured 
or  calculated  waveforms  or  both. 

As  a  demonstration  of  capability,  the  computer  program  was  used  to 
statistically  analyze  selected  data  from  the  DISC  Test  I  event  and  objectively 
compare  particle  velocity  measurements  made  in  DISC  Test  II  with  the  expected 
value  waveform  obtained  from  probabilistic  prediction  calculations. 

5 . 2  RECOMMENDATIONS 

It  is  recommended  that  the  objective  discrepancy  measures  examined  in 
this  report  be  used  whenever  comparisons  of  two  or  more  waveforms  are  made. 

It  is  also  recommended  that  the  technique  be  extended  to  objectively  quantify 
differences  in  laboratory-  and  field-generated  material  property  test  results. 
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APPENDIX  A 


USER'S  GUIDE  FOR  COMPUTER  PROGRAM  WCT 


A. 1  INTRODUCTION 

This  user's  guide  for  the  computer  program  WCT  describes  typical  input 
and  output,  contains  a  glossary  of  the  variables,  a  flow  chart,  and  a  listing 
of  the  program,  and  presents  sample  tabulated  output  from  two  example  runs. 
Program  WCT  has  been  coded  in  Honeywell  Level  66  Fortran  for  the  timesharing 
subsystem  of  the  Honeywell  DPS-1  digital  computer  currently  operated  at  WES. 

A. 2  INPUT 

The  digitized  waveforms  for  input  to  program  WCT  can  be  read  directly 
from  tapes  or  can  be  written  to  binary  data  files  for  subsequent  access 
through  the  DPS-1  timesharing  subsystem.  Each  input  waveform  consists  of  two 
records  which  are  treated  as  one  tape  file  by  WCT.  The  first  record  contains 
the  number  of  digitized  data  points,  the  digital  time  increment,  and  an  iden¬ 
tification  (ID)  label.  For  measured  waveforms,  this  ID  label  consists  of 
three  20-character  alpha-numeric  variables  that  contain  all  pertinent  informa¬ 
tion  about  the  gage  and  type  of  data;  for  calculated  waveforms,  the  ID  label 
consists  of  a  title  up  to  60  characters  in  length.  The  record  contains  the 
digitized  waveform  as  a  single  array  of  data  points.  All  other  input  is  in 
free-field  format,  and  the  program  will  call  for  the  variables  by  name. 

The  first  line  of  input  variables  contains  the  following  information: 

XFINAL - Final  time  for  calculations,  msec. 

DX - - — Time  increment,  msec. 

NTPLOT - Type  of  calculations  desired:  For  NTPLOT  *  1,  objective  discrepancy 

measures  (Equations  2.5  through  2.14  of  Chapter  2)  are  computed; 


for  NTPLOT  =  2,  expected  value  and  standard  deviation  are  computed; 
and  for  NTPLOT  =  3,  expected  value,  standard  deviation,  and 
derivatives  with  respect  to  time  for  expected  value  and  standard 
deviation  are  computed. 

SEARV - Search  value  for  normalizing  the  arrival  times  of  the  records 

(about  1/2  percent  of  peak  value  of  the  data) .  If  SEARV  is  input 
as  zero,  the  data  will  be  read  from  the  beginning  of  the  record. 

ISKIP - Print  SKIP  increment. 

NIBASE - Number  of  integrations  for  base  record. 

NICOMP - Number  of  integrations  for  comparison  records.  (The  base  record 

is  treated  separately  from  the  other  records.  This  will  allow  a 
comparison  of  the  base  record  to  one,  or  more,  measured  or  calcu¬ 
lated  waveforms . ) 

The  second  line  of  input  contains  the  following  variables  which  are  required 

by  the  program  for  every  record: 

NSORCE -  =  0;  no  more  waveforms  to  be  read  in. 

=  1;  measured  waveforms. 
m  2;  calculated  waveforms. 

NFILE - File  number.  If  NFILE  *  0,  the  program  will  ask  for  the  name  of 

a  data  file  containing  the  waveform(s). 

The  third  line  of  input  is  the  name  of  the  data  file,  called  "FILE,"  contain¬ 
ing  the  waveform.  Finally,  the  program  asks  for  the  value  of  ANS,  which  gives 

the  user  options  to  obtain  plots  or  tables  or  both,  as  described  below. 


A.  3  OUTPUT 


WCT  output  consists  of  optional  time  history  plots  or  tabulated  data  or 
both.  In  addition,  the  input  records  can  be  plotted  either  before  or  after 
preprocessing  or  both.  A  table  of  maximum  and  minimum  values  if  produced  for 
all  computations. 

The  type  of  output  depends  on  the  value  of  NTPLOT.  For  NTPLOT  =  1,  the 
output  consists  of  the  magnitude  error  (Equation  2.8),  the  phase-and-f requency 
error  (Equation  2.9),  and  the  combined  error  (Equation  2.10).  If  there  is 
more  than  one  comparison  record,  the  mean  of  each  error  (i.e.,  ensemble 
averages.  Equations  2.12,  2.13,  and  2.14)  is  also  computed.  For  NTPLOT  =  2, 
the  output  consists  of  expected  (or  mean)  values  and  standard  deviation.  If 
the  arrival  times  of  the  records  have  been  normalized,  the  expected  waveforms 
can  be  plotted  against  the  time  associated  with  expected  arrival  time  and  the 
expected  waveform  plus  or  minus  standard  deviation  can  be  plotted  against  the 
expected  arrival  time  plus  or  minus  the  standard  deviation,  respectively. 

For  NTPLOT  *»  3,  the  output  consists  of  data  identical  to  NTPLOT  =  2,  plus  the 
derivatives  with  respect  to  time  of  (a)  the  expected  value,  (b)  the  expected 
value  plus  one  standard  deviation,  and  (c)  the  expected  value  minus  one 
standard  deviation. 

A. 4  GLOSSARY 
A. 4.1  Main  Program 

ANS  Character  variable  through  which  the  types  of  outputs  are 

chosen. 

CEF 

DE(I) 


Combined  error  factor  (Equation  2.7). 

Derivative  with  respect  to  time  of  E(I)  at  the  Ith  time 
step. 

A3 


DEM(I) 

DEMAX 

DEMMIN 

DEP(I) 

DEPMAX 

DT(J) 

DX 

DXI 

DX02 

E(I) 

ECMN(K) 

ECMX(K) 

ECOM(I.K) 

ECOMAV(I) 

EM(I) 


Derivative  with  respect  to  time  of  EM(I)  at  the  Ith  time 
step. 

Maximum  value  of  DE(I). 

Minimum  value  of  DEM(I). 

Derivative  with  respect  to  time  of  EP(I)  at  the  Ith  time 
step. 

Maximum  value  of  DEP(I). 

Time  increment  for  the  Jth  record,  msec. 

Time  increment  for  calculations,  msec. 

1/DX. 

DX/2. 

Expected  value  (mean  of  given  set  of  records  at  the  Ith 
time  step. 

Minimum  value  of  ECOM(I,K)  at  the  Kth  record. 

Maximum  value  of  ECOM(I,K)  at  the  Kth  record. 

Combined  error  between  the  Kth  record  and  the  base  one  at 
the  Ith  time  step  (Equation  2.10). 

Mean  combined  error  at  the  Ith  time  step  (Equation  2.14). 
E(I)  Minus  one  standard  deviation  at  the  Ith  time  step. 


EMAG(I,K) 


Magnitude  error  between  the  Kth  record  and  the  base  one 
at  the  Ith  time  step  (Equation  2.8). 


EMAGAV(I) 


Mean  magnitude  error  at  the  Ith  time  step  (Equation  2.12). 


EMAX 

EMIN 

EMMIN 

EMMN(K) 

EMMX(K) 

EP(I) 

Er  HS  ( I ,  K) 

EPHSAV(I) 

EPMAX 

EPMN(K) 

EPMX(K) 

ES 

11(1) 

ICM1 

ICNT 

ISKIP 


Maximum  value  of  E(I). 

Minimum  value  of  E(I). 

Minimum  value  of  EM(I). 

Minimum  value  of  EMAG(I,K)  at  the  Kth  record. 

Maximum  value  of  EMAG(I,K)  at  the  Kth  record. 

E(I)  plus  one  standard  deviation  at  the  Ith  time  step. 

Phase-and-frequency  error  between  the  Kth  record  and  the 
base  one  at  the  Ith  time  step  (Equation  2.9). 

Mean  phase-and-frequency  error  at  the  Ith  time  step 
(Equation  2.13). 

Maximum  value  of  EP(I). 

Minimum  value  of  EPHS(I,K)  at  the  Kth  record. 

Maximum  value  of  EPHS(I,K)  at  the  Kth  record. 

Temporary  variable  for  computing  the  expected  value. 
Integral  of  the  base  record  squared  at  the  Ith  time  step. 
Number  of  records  to  be  compared  with  the  base  record. 
Total  number  of  records  to  be  processed. 

Print  SKIP  increment. 


A5 


MAXC(K) 

MAXM 

MCF 

N1 

NC 

NIBASE 

NICOMP 

NINT 

NINVERS 

NP 

NPOINT 

NPIS(J) 

NTPLOT 

PCF 

PEF(K) 


One  plus  maximum  absolute  value  of  the  record  at  K. 

Maximum  absolute  of  the  base  record. 

Magnitude  correlation  factor  (Equation  2.5). 

Parameter  variable  for  setting  the  maximum  number  of 
comparison  cases  to  be  processed. 

Parameter  variable  for  setting  the  maximum  number  of 
records  to  be  processed. 

Number  of  times  for  which  the  base  record  must  be 
integrated. 

Number  of  times  for  which  the  waveforms  must  be  integrated. 
Temporary  counter  for  NIBASE  and  NICOMP. 

1/lCMl. 

Parameter  variable  for  setting  the  maximum  number  of  time 
steps  that  can  be  processed. 

Number  of  time  steps  to  be  used  for  given  calculations. 
Number  of  time  steps  in  the  Jth  record. 

Variable  determines  the  desired  type  of  output. 
Phase-and-frequency  correlation  factor  (Equation  2.6). 

Peak  error  factor  of  the  Kth  record. 


PL0T2 


A  subroutine  for  plotting  on  a  Tektronix  4662  interactive 
digital  plotter  (Note:  It  is  not  the  intent  of  this  user’s 
guide  to  explain  the  use  of  PL0T2) . 


RICNT 


l./ICNT. 


RNM1 

SS 

ST 

SUM1 

SUM2 

TAR(K) 

TE(I) 

TE1 

TITLE (K) 

TM(I) 

TM1 

TP(I) 

TP1 

X(J,K) 

XCUR 

XFINAL 


1/ (ICNT-1) . 

Variance. 

Standard  deviation. 

Temporary  variable  for  computing  MCF. 

Temporary  variable  for  computing  PCF. 

Arrival  time  for  the  Kth  record,  msec. 

Time  associated  with  E(I)  at  the  Ith  time  step,  msec. 

Expected  arrival  time  of  the  given  set  of  records,  msec. 

Title  with  up  to  60  characters  for  identifying  the  Kth 
record. 

Time  associated  with  EM(I)  at  the  Ith  time  step,  msec. 
TE1-ST. 

Time  associated  with  EP(I)  at  the  Ith  time  step,  msec. 

TE1+ST. 

Time  associated  with  the  Kth  record  at  the  Jth  time  step, 
msec. 

Temporary  variable  used  for  interpolation. 

The  length  of  the  calculation,  msec. 


XX(I) 


YNEXT.YT 


Time  at  the  Ith  time  step,  msec. 


An  array  for  storing  the  input  data. 


Temporary  variables  for  integration. 


An  array  for  storing  the  preprocessed  data. 


A. 4. 2  Subrouting  READ IN 


ATTACH 


BCDASC 


C1.C2.C3 


DETACH 


DT(ICNT) 


ISKIP 


Nl.NC.NP 


NFILE 


System  subroutine  for  opening  a  permanent  file. 

System  subroutine  for  converting  from  BCD  to  ASCII. 

Identification  for  measured  data;  BCD  labels  (converted 
to  ASCII  for  title). 

System  subroutine  for  closing  a  permanent  file. 

Time  increment  (in  seconds,  converted  to  milliseconds) 
for  the  ICNT  record. 

Time  increment,  msec. 

Variable  name  for  input  file. 

Record  counter. 

Print  SKIP  increment. 

(See  main  program) . 

Tape  file  number  in  the  input  file. 


NIBASE,  NICOMP  (See  main  program). 


NPOINT 


NPS 

NPT 

NPTX(ICNT) 

NSORCE 

NSTRT 

NTPLOT 

I 

SEARV 

i 

i 

TAR(ICNT) 

TDUM 

TITLE (ICNT) 


Number  of  points  to  be  used  for  calculations.  NPOINT  will 
be  reduced  if  any  input  record  contains  fewer  than  NPOINT 
data  points  (this  would  also  reduce  the  value  of  XFINAL) . 

Maximum  number  of  data  points  from  an  input  record  to  be 
searched  in  order  to  obtain  the  arrival  time. 

Total  number  of  data  points  to  be  read  from  an  input  record 

(See  main  program.) 

Origin  of  input  data:  NSORCE  =  1  for  measured  data; 

NSORCE  =  2  for  calculated  data. 

The  number  of  time  steps  at  which  the  signal  arrives. 

(See  main  program.) 

(See  the  input.) 

Arrival  time  (in  seconds,  converted  to  milliseconds)  for 
record  number  ICNT. 

Identification  label.  BCD  label  (converted  to  ASCII  for 
title) . 

(See  main  program.) 


X,XFINAL,XX,Y,YY  (See  main  program.) 


30  NINT  *  NICOMP 


INTERPOLATE  Y(I,  K)  INTO  YY(I,  K) 


CALCULATE  EMAG(I,  K); 
EPHS(I,  K);  tCOM(I,  K) 


T 


F 


CALCULATE  EMAGAV(I); 
EPHSAV(I);  ECOMAV(I) 


PRINT:  WANT  TABULATED  OUTPUT? 
READ  ANS 


TABULATE  THE  MAXIMUM  AND  MINIMUM  VALUES 


TABULATE  THE  CALCULATED  DATA 


PRINT:  WANT  TABULATED  OUTPUT? 


TABULATE  THE  MAXIMUM  AND  MINIMUM  VALUES 


TABULATE  THE  CALCULATED  DATA 

* 

PRINT:  WANT  PLOTS? 

READ  ANS 


PLOT  THE  CALCULATED  DATA 


PLOT  THE  CALCULATED  DATA 


PRINT:  READ  XFINAL,  DX,  NTPLOT,  SEARV, 
ISKIP.,  NIBASE,  NICOMP 

READ  XFINAL,..., NICOMP 


DECREMENT  ICNT  BY  1 
COMPUTE  TIME  ARRAYS 


v.7 


A.  6  EXAMPLES 


OLD  WCT 
IFRN 


12/08/82  21.070 


READ  XFINAL»DX»NTPLOT » SEARV* ISKIP.NIBASEtNICOHP 

=7  ,02  1  .5  50  0  0 

NPOINT  =  351 

READ  NSORCErNFILE 

=1  0 

INPUT  FILE  ? 

*R0SD222/0»JTB1 
READ  NFILE 

=i 

READ  NSORCErNFILE 

=1  2 

READ  NSORCErNFILE 
=1  5 

READ  NSORCErNFILE 
=1  3 

READ  NSORCEf NFILE 
=0  0 

WANT  PLOT  OF  INPUT  DATA? 

®N0 

WANT  PLOT  OF  PREPROCESSED  DATA? 

=N0 

WANT  TABULATED  OUTPUT  ? 

=YES 


CASE 

1 

2 

NAXC  = 

41.598 

34.930 

rt'F  = 

-0.097 

-0.241 

ENHX  * 

105.420 

88.186 

ENMN  * 

-0.460 

-0,043 

EPHX  = 

0.987 

0.995 

EPNN  * 

0. 

0. 

ECMX  ■ 

105.423 

88.191 

ECHN  = 

-0.777 

-0.730 

BASE 

3.5-0-0-AB 

2-6 

CASE 

1  6-0-45-A8 

2-7 

I 

EHAG 

EPHS 

2 

-0.460 

0.626 

50 

105.243 

0.898 

100 

0.148 

0.777 

150 

0.093 

0.647 

200 

0.097 

0.588 

250 

0.069 

0.537 

300 

0.073 

0.497 

350 

0.059 

0.470 

CASE 

2  2-0-300-AB 

2-l< 

3 

19.221 

-0,583 

51.388 

-0.321 

0.991 

0. 

51.398 

-0.704 


PRESSURE  KPA  DISC  TEST  1 

PRESSURE  KPA  DISC  TEST  1 

ECON 

-0.777 

105.247 

0.791 

0.654 

0.596 

0.542 

0.503 

0.474 

PRESSURE  KPA  DISC  TEST  1 


DISC  TEST  1 


50 

85.993 

0.931 

85.998 

100 

0.000 

0.733 

0.733 

150 

-0.012 

0.583 

-0.583 

200 

-0.007 

0.530 

-0.530 

250 

-0.033 

0.486 

-0.487 

300 

-0.031 

0.450 

-0.451 

350 

-0.043 

0.425 

-0.427 

CASE 

3  4-0-90-AB 

2-8 

PRESSURE 

I 

EMAG 

EPHS 

ECOH 

2 

3.910 

0.592 

3.955 

50 

43.846 

0.876 

43.854 

100 

-0.293 

0.503 

-0.582 

150 

-0.307 

0.411 

-0.513 

200 

-0.238 

0.366 

-0.437 

250 

-0.248 

0.330 

-0.413 

300 

-0.223 

0.302 

-0.375 

350 

-0.219 

0.279 

-0.355 

I 

EMAGAV 

EPHSAV 

ECOMAV 

2 

1.347 

0.598 

1.852 

50 

78.360 

0.902 

78.366 

100 

-0.048 

0.671 

-0.702 

150 

-0.075 

0.547 

-0.583 

200 

-0.049 

0.494 

-0.521 

250 

-0.071 

0.451 

-0.481 

300 

-0.060 

0.416 

-0.443 

350 

-0.068 

0.391 

-0.419 

5 


WANT  PLOTS  ? 

=NO 

PTU-SEC  *  1.85 


12/08/82 


21.100 


READ  XFINAL » DX  t NTPLOT r  SEARV » I SKIP. NI BASE » NI COMP 
=7  .02  3  ,5  50  1  1 
NPOINT  =  351 

READ  NSORCE,NFILE 
=1  0 

INPUT  FILE  ? 

=R0SD222/QUTB291 
READ  NFILE 

=i 

READ  NSORCE'NFIIE 

=1  2 

READ  NSORCErNFILE 
=1  5 

READ  NSORCEf NFILE 
=1  3 

READ  NSORCE'NFILE 
=0  0 

WANT  PLOT  OF  INPUT  DATA? 

=N0 

WANT  PLOT  OF  PREPROCESSED  DATA? 

=N0 

WANT  TABULATED  OUTPUT  ? 

=YES 


ENX 

EHN 

EPNX 

EHHN 

41.336 

0. 

44.986 

-0.008 

TE1 

TP1 

TH1 

0.960 

1.646 

0.274 

CASE 

1 

2 

3 

TAR  = 

0. 

1.600 

1.000 

I 

E 

EP 

EH 

50 

7.048 

12,494 

1.602 

100 

14.246 

18.197 

10.295 

150 

20.785 

24.927 

16.643 

200 

26.490 

30.624 

22.356 

250 

31.477 

35.307 

27.647 

300 

36.760 

40.705 

32.814 

350 

41.255 

44.908 

37.601 

WANT  PLOTS  ? 

*N0 

WANT  TABULATED  OUTPUT  ? 
BYES 


DENX 

18.408 


DEMN 

0. 


DEPHX 

31.840 


DENNN 

-0.880 


A. 7  PROGRAM  LISTING 


1000*#RUNH*IR0SD441/PL0TS*R 


1010C 

1020C 

1030C 

1040C 

1050C 

1060 

1070 

1080 

1090 

1100 

1110 

1120 

1130 

1140 

1150 

1160 

1170 

1180 

1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260C 

1270 

1280C 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400C 

1410C 

1420C 

1430C 

1440 

1450 

1460 

1470 

1480 

1490 

1500 


PROGRAM  WCT 

CALCULATIONS  OF  STATISTICAL  MEASURES  FOR 
COMPARISON  OF  WAVEFORMS 

PARAMETER  NC  *  10»N1  *  NC-1»NP  =  200 
REAL  IlfMCF'HAXC'MAXHfNINVRS 
CHARACTER  TITLE  *60»ANS*1 

DIMENSION  II (NP) »EMAGAV(NP) »EPHSAV(NP) »ECOMAV(NP) t 
t  EMAG(NPfNl ) »EPHS(NP»N1) »EC0M(NP»N1 ) t 

l  TE(NP) >TP(NP) »TM(NP) »E(NP) tEP(NP) rEM(NP) » 

t  DE(NP) »DEP(NP>  fDEM(NP) 

COMMON  /INPUT/  XFINAL»DX»NTPLOT»ISKIP»NIBASE»NICOMP, 
t  ICNT»NPOINT » DT(NC) »TAR(NC) f NPTS(NC) 

COMMON  /ARRA1/  X(NP»NC) »Y<NP»NC> rXX(NP) »YY<NP»NC> » 

S  MAXC ( N 1 >  *  PEF  <  N 1 ) t EMMX ( N 1 ) r EMMN ( N1 ) i EPMX ( N 1 ) r 

S  EPNN(N1 ) »ECMX(N1 ) »ECMN(N1) > TITLE (NC) 

EQUIVALENCE  (EMAGAVd > » II ( 1 ) »X< 1 r NC) ) » (EMAGCl 1 1 ) »X( 1 1 1 ) >  t 
t  (EPHSAVI l>#Y(lf NC>)» (EPHS(1 »1 ) » Y(1 r 1 ) >  t 

t  (ECOMAVd)  rYY(lfNC) )  t  (ECOMd f  1 ) rYY(l r  1  > ) 

EQUIVALENCE  (TEd )  »Xd » 1 ) )  t  (TPd  >  »Xd  ,2)  >  f  (TMd)  f  Xd  »3)  >  * 

S  (Ed)fYdvl))»(EPd)fYd»2))i(EN<l)iYd»3))f 

t  (DE(l)fYYdrl)) » (DEPd )  i  YYd  »2) )  r  (DEMd)  fYY(l?3>) 

CALL  PTIME(PTI) 

CALL  FPARAM(1»80) 

CALL  READIN 

PRINT f 'WANT  PLOT  OF  INPUT  DATA?* 

READfANS 

IF( ANS.NE.1HY)  60  TO  10 
1  CONTINUE 
DO  5  K  =1 » ICNT 

CALL  PL0T2(X(1»K) »Y(1»K>»NPTS(K) ) 

5  CONTINUE 

PRINT» 'WANT  REPLOT  ?• 

READfANS 

IF( ANS.EQ. 1HY)  GO  TO  1 
10  CONTINUE 

PERFORM  INTEGRATIONS  AS  NEEDED  ON  DATA  TO  OBTAIN 
DESIRED  QUANITIES  FOR  COMPARISON 

NINT  *  NIBASE 
DO  30  J  =1 t ICNT 
IF(NINT.EQ.O)  60  TO  30 
DO  20  L  *1 fNINT 
DX02  =  .5*DT(J) 

YNEXT  *  Yd » J>  +  Y(2»J) 

Y(1fJ)  =  0. 


1510  Y (NPTS<  J)  +  1  * J)  =  0. 

1520  DO  20  I  =2*NPTS< J) 

1530  YT  *  YNEXT 

1540  YNEXT  »  Y<I»J)  +  YU  +  1*J> 

1550  Y < I » J )  =  Y(I— 1#J)  +  DX02  *  YT 

1560  20  CONTINUE 

1570  30  NINT  =  NICONP 

1580C 

1590C  INTERPOLATE  VERTICAL  ARRAYS  FOR  SAKE  DX  IF  REQUIRED 

1600C 

1610  DO  60  K  =1 *  ICNT 

1620  IF(DT(K) «LE«  DX)  60  TO  55 

1630  YY( 1 »K)  *  Y(1*K) 

1640  XCUR  -  DX 

1650  I  =  1 

1660  DO  50  J  =2*NPTS(K) 

1670  40  IF(X(J*K)»LT« XCUR )  GO  TO  50 

1680  1=1+1 
1690  JH1  3  J  -  1 

1700  YY< I »K)  =  Y( JH1 *K)  +  <Y< J*K)-Y< JN1 *K) )*<XCUR-X( JM1 *K> )/ 

1710  i  (X( J*K)-X( JM1*K) ) 

1720  XCUR  *  I  t  DX 

1730  IF( I-NPOINT)  40»60»60 

1740  50  CONTINUE 

1750  GO  TO  60 

1760  55  CONTINUE 

1770C 

1780C  IF  INTERPOLATION  NOT  REQUIRED  TRANSFER  Y  ARRAY  INTO  YY  ARRAY 
1790C 

1800  DO  58  I  =1 »NPOINT 
1810  58  YY( I »K)  =  Y(I»K) 

1820  60  CONTINUE 

1830C 

1840C  SET  UP  HORIZONTAL  ARRAY 

1850C 

1860  DO  70  I  =1 »NPOINT 

1870  XX<I>  *  DX  *  <1-1 » 

1880  70  CONTINUE 

1890C 

1900  PRINT* 'WANT  PLOT  OF  PREPROCESSED  DATA?* 

1910  READrANS 

1920  IF(ANS.NE.IHY)  GO  TO  90 

1930  80  CONTINUE 

1940  DO  85  K  =1*ICNT 

1950  CALL  PL0T2(XX* YY ( 1 *K) iNPOINT) 

1960  85  CONTINUE 

1970  PRINT* ‘WANT  REPLOT  ?■ 

1980  READrANS 

1990  IF(ANS.EQ.IHY)  GO  TO  80 

2000  90  CONTINUE 

2010  IF(NTPLOT.GT.l)  GO  TO  400 

2020C 

2030C  FORM  INTEGRALS  FOR  CORRELATIONS 


2050 

2060 

2070 

2080 

2090 

2100 

2110 

2120 

2130 

2140 

2150 

2160 

2170 

2180 

2190 

2200 

2210 

2220 

2230 

2240 

2250 

2260 

2270 

2280 

2290 

2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2450 

2460 

2470 

2480 

2490 

2500 

2510 

2520 

2530 

2540C 

2550C 

2560C 


11X02  =  0.54DX 
ICM1  =  ICNT-1 
11(1)  =  0. 

MAXM  =  ABS(YYdd)) 

DO  100  I  =2»NP0INT 

MAXM  =  MAX( A8S( YY ( I r 1 ) > rNAXM) 

11(1)  =  11(1-1)  +  DX02  *  ( YY< I » 1 >**2+YY( 1-1 » 1 >**2> 

100  CONTINUE 
K  =  2 

110  KH1  =  K  -  1 
MCF  =  1.0 
PCF  =  1.0 
EHHX(KHl)  =  0. 

ENNN(KHl)  =  0. 

EPHX(KMl)  =  0. 

EPNN(KHl)  =  0. 

MAXC(KMl)  =  ABS(YY  < 1 »K> ) 

PEF(KMl)  =  0. 

EMAG(1»KM1)  =  0. 

EPHS(IfKMI)  =  0. 

SUN1  =  0. 

SUM2  =  0. 

DO  130  I  =2»NP0INT 

NAXC(KHl)  =  MAX(MAXC(KM1 ) » ABS(YY  < I » K ) ) ) 

SUN1  *  SUm  +  DX02  *  (YY(I-l»K>M2m(IfK)«2) 

SUM2  =  SUN2  +  DX02  *  < YY( 1-1 d )*YY( 1-1 »K)+YY( I t 1 )*YY( I »K) ) 

MCF  =  SORT ( MAX (SUM1 » .001 )/MAX( II ( I) » .001 ) ) 

PCF  a  MAX(ABS(SUM2) > .001 )  /  MAX(SQRT(SUM1*I1 ( I ) ) » .001 ) 
EMAG(IfKMl)  =  (MCF-i . ) 

EPHS(IfKMl)  =  (l.-PCF) 

EMMX(KMl)  =  MAX(EMMX(KM1 ) r EMAG( I»KM1 ) ) 

EPMX(KMl)  =  MAX(EPMX(KM1)»EPHS(IfKM1)) 

EMMN(KMl)  »  MIN(EMMN(KM1 ) »ENAG( I »KMl ) ) 

EPMN(KMl)  =  MIN(EPMN(KM1 ) »EPHS( I»KM1 ) > 

130  CONTINUE 

PEF(KMl)  =  MAXC(KM1)/MAX(MAXMf.001)-1. 

K  “  K  ♦  1 

IF(K.LE.ICNT)  GO  TO  110 
DO  140  K  =  1  dCMl 
ECMX(K)  =  0. 

ECMN(K)  =  0. 

ECOMd  »K)  =  0. 

DO  140  I  »2fNP0INT 

CEF  *  SQRT((EMAG(I»K))»*2HEPHS(IfK))**2) 

ECOM(I>K)  -  SIGN(CEF»EMAG( I»K) ) 

ECMX(K)  =  MAX(ECMX(K)rECOM(I»IO) 

ECMN(K)  «  MIN(ECMN(K)»ECOH(IfK)) 

140  CONTINUE 

IF(ICNT.EQ.2)  GO  TO  200 

IF  MORE  THAN  1  COMPARISON  COMPUTE  THE  AVERAGES 


VIV.VJS.* 


2570  EMAGAVd)  =  0. 

2580  EPHSAVd)  *  0. 

2590  ECOHAVd)  =  0. 

2600  NINVRS  =  1./ICH1 

2610  DO  160  I  a2»NP0INT 

2620  EHAGAVd >  =  0. 

2630  EPHSAVd)  a  0. 

2640  ECOHAVd)  =  0. 

2650  DO  150  K  al.JCHl 

2660  EHAGAVd)  =  EHAGAVd)  +  EHAGd.K) 

2670  EPHSAVd)  a  EPHSAVd)  +  EPHSd.K) 

2680  ECOHAVd)  a  ECOHAVd)  +  ABS(ECOHd.K) ) 

2690  150  CONTINUE 

2700  EHAGAVd)  *  EHAGAVd)  *  NINVRS 

2710  EPHSAVd)  =  EPHSAVd)  *  NINVRS 

2720  ECOHAVd)  »  SIGN ( NINVRStECOHAVd )  .EHAGAVd ) ) 

2730  160  CONTINUE 

2740C 

2750  200  CONTINUE 

2760C  OUTPUT  PHASE 

2770C 

2780  PRINT. 'WANT  TABULATED  OUTPUT  ?* 

2790  READ. ANS 

2800  K1  *  1 

2810  K2  =  ( I CM  1/8 )  ♦  1 

2820  K3  =  HINI7.ICH1) 

2830  DO  205  I  al.K2 

2840  PRINT  665. (K*KaKl»K3) 

2850  PRINT  670. (HAXC(K) .K-Kl .K3) 

2860  PRINT  671 . (PEF(K) »KaKl »K3) 

2870  PRINT  672. (EHHX(K) »K=K1 »K3> 

2880  PRINT  673. (EHHN<K) »KaKl »K3) 

2890  PRINT  674. (EPMX(K) . K=K1 .K3) 

2900  PRINT  675. (EPHNIK) »KaKl »K3) 

2910  PRINT  676.<ECHX(K)»KaKl.K3) 

2920  PRINT  677. (ECHN(K) .KaKl .K3) 

2930  K1  =  8 

2940  K3  =  ICH1 

2950  205  CONTINUE 

2960  IF(ANS.NEdHY)  GO  TO  300 

2970  PRINT  650. TITLE < 1 ) 

2980  DO  220  Kal.ICHl 

2990  PRINT  660.K.TITLE(K+l) 

3000  PRINT  600 

3010  PRINT  610.2.EHAG(2.K) »EPHS(2.K) »EC0H(2»K) 

3020  DO  210  I  =ISKIP.NPOINT . ISKIP 

3030  PRINT  610. 1 .EHAGd  . K ) . EPHSd  .K)  .ECOHd  .K) 

3040  210  CONTINUE 

3050  220  PRINT. 

3060  IF(ICNT.EQ.2)  GO  TO  260 

3070  PRINT  620 

3080  PRINT  630. 2. EHAGAVd). EPHSAV(2).EC0HAV<2> 

3090  DO  230  I  *ISKIP.NPOINT. ISKIP 
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3100  PRINT  630 t I t EHAGAV ( I ) tEPHSAV( I >  t EC0MAV( I ) 

3110  230  CONTINUE 

3120  PRINT  t 

3130  260  CONTINUE 

3140  PRINT t 'WANT  PLOTS  ?• 

3150  READ* ANS 

3160  IF(ANS.NE.IHY)  60  TO  99? 

3170  300  CONTINUE 

3180  PRINT t 'WANT  PLOT  OF  XX-EHAG  ?' 

3190  READtANS 

3200  IF(ANS.NE.IHY)  GO  TO  340 

3210  DO  335  K  =ltICHl 

3220  CALL  PL0T2(XXtENAG< 1 tK) tNPQINT) 

3230  335  CONTINUE 

3240  340  PRINT t ‘WANT  PLOT  OF  XX-EPHS  ?* 

3250  READtANS 

3260  IF(ANS.NE.IHY)  GO  TO  350 

3270  DO  345  K  =ltICNl 

3280  CALL  PL0T2(XXt EPHS( 1 »K> tNPQINT) 

3290  345  CONTINUE 

3300  350  PRINT » 'WANT  PLOT  OF  XX-ECOH  ?' 

3310  READtANS 

3320  IF(ANS.NE.IHY)  GO  TO  360 

3330  DO  355  K  =ltICNl 

3340  CALL  PL0T2(XXtEC0M( 1 tK) t NPOINT) 

3350  355  CONTINUE 

3360  360  IF < ICNT . EQ« 2 )  GO  TO  380 

3370  PRINT t ‘WANT  PLOT  OF  XX-EHAGAV  ?• 

3380  READtANS 

3390  IF(ANS.EQ.IHY)  CALL  PL0T2 ( XX t EHAGAV t NPOINT) 

3400  PRINT t 'WANT  PLOT  OF  XX-EPHSAV  ?' 

3410  READtANS 

3420  IF(ANS.EQ.IHY)  CALL  PL0T2<XXtEPHSAVtNP0INT) 

3430  PRINT t ‘WANT  PLOT  OF  XX-ECOHAV  ?' 

3440  READtANS 

3450  IF(ANS.EQ.IHY)  CALL  PL0T2<XXt ECOMAVtNPOINT) 

3460  380  PRINT t  'WANT  REPLOT  ?* 

3470  READtANS 

3480  IF(ANS.EQ.IHY)  GO  TO  300 

3490  GO  TO  999 

3500C 

3510C  CALCULATE  NEAN  AND  STANDARD  DEVIATIONS 

3520C 

3530  400  CONTINUE 

3540  RICNT  =  l./ICNT 

3550  RNN1  =  1 ./( ICNT-1 ) 

3560  ES  =  0. 

3570  EMAX  =  0. 

3580  ENIN  =  0. 

3590  EPNAX  =  0. 

3600  ENHIN  =  0, 

3610  DO  410  K*1 t ICNT 

3620  ES  8  ES  f  TAR(K) 


3630  410  CONTINUE 

3640  TE1  =  ESfRICNT 

3650  SS  =  0. 

3660  DO  420  K=1»ICNT 

3670  SS  =  SS  4  <TAR(K)-TE1>#*2 

3680  420  CONTINUE 

3690  ST  =  SORT (SStRNMl ) 

3700  TP1  =  TElfST 

3710  TNI  =  TE1-ST 

3720  DO  450  I=1»NP0INT 

3730  ES  =  0. 

3740  DO  430  K*1»ICNT 

3750  ES  =  ES  4  YY(I»K) 

3760  430  CONTINUE 

3770  E(I>  =  EStRICNT 

3780  SS  *  0. 

3790  DO  440  K=1»ICNT 

3800  SS  =  SS  4  <YY(I.K)-E(I))**2 

3810  440  CONTINUE 

3820  ST  =  SORT (SSIRNN1 > 

3830  £P(I )  =  E(I)4ST 

3840  EM ( I )  =  E  C I > —ST 

3850  EMAX  =  MAX ( E ( I ) »ENAX) 

3860  EMIN  =  MIN(E( I ) >EMIN) 

3870  EPMAX  =  MAX<EP( I) »EPMAX) 

3880  EMMIN  *  MINIEHI I ) ,EMHIN> 

3890  TE(I)  =  XX(I)4TE1 

3900  TP ( I >  =  XX< I )4TP1 

3910  TM< I >  *  XX<I)4TM1 

3920  450  CONTINUE 

3930  PRINT » 'WANT  TABULATED  OUTPUT  ?' 

3940  READ»ANS 

3950  PRINT  680>EMAXf ENIN 'EPMAX 'EMMIN 'TElrTPl r TNI 

3960  K1  =  1 

3970  K2  =  ( ICNT/8)  4  1 

3980  K3  *  MIN(7' ICNT ) 

3990  DO  455  I  *1»K2 

4000  PRINT  665r (K'KSK1'K3) 

4010  PRINT  682' (TAR(K) 'KSK1 'K3> 

4020  K1  =  8 

4030  K3  *  ICNT 

4040  455  CONTINUE 
4050  PRINT' 

4060  IF(ANS.NE.IHY)  GO  TO  470 

4070  PRINT  690 

4080  DO  460  I  «ISKIP'NPOINT» ISKIP 

4090  PRINT  630'I»E<I)'EP(I)'EM(I) 

4100  460  CONTINUE 
4110  PRINT' 

4120  PRINT' 

4130  PRINT » ‘WANT  PLOTS  ?• 

4140  READ'ANS 

4150  IF(ANS.NE.IHY)  GO  TO  500 
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4160  470  CONTINUE 

4170  PRINT r 'WANT  PLOT  OF  TE-E  ?* 

4180  READrANS 

4190  IF( ANS.EQ. 1HY)  CALL  PL0T2<TErErNP0INT> 

4200  PRINT r ‘WANT  PLOT  OF  TP-EP  ?* 

4210  READrANS 

4220  IFtANS.EQ.lHY)  CALL  PL0T2<TPr EPrNPOINT) 

4230  PRINT r ’WANT  PLOT  OF  TH-EN  ?* 

4240  READiANS 

4250  IF(ANS.EQ.IHY)  CALL  PL0T2<TNrENrNPDINT) 

4260  PRINT r 'WANT  PLOT  OF  XX-E  ?* 

4270  READrANS 

4280  IF(ANS.EO.IHY)  CALL  PL0T2(XXtErNP0INT) 

4290  PRINT » 'WANT  PLOT  OF  XX-EP  ?• 

4300  READrANS 

4310  IF(ANS.EQ* 1HY>  CALL  PL0T2(XXrEPrNP0INT) 

4320  PRINT r ‘WANT  PLOT  OF  XX-EM  ?* 

4330  READrANS 

4340  IF(ANS.EQ.IHY)  CALL  PL0T2(XXrEHrNP0INT> 

4350  PRINT r ’WANT  REPLOT  ?* 

4360  READrANS 

4370  IF(ANS.EQ.IHY)  GO  TO  470 

4380  500  CONTINUE 

4390  IF(NTPLOT . LT*3>  GO  TO  999 

4400C 

4410C  CALCULATE  DERIVATIVE  WITH  RESPECT  TO  TINE 

4420C  FOR  NEAN  AND  STANDARD  DEVIATIONS 

4430C 

4440  DENAX  =  0. 

4450  DENIN  =  0. 

4460  DEPNAX  *  0. 

4470  DENMIN  =  0. 

4480  DXI  *  l./DX 

4490  DE( 1 )  =  0. 

4500  DEP(l)  =  0. 

4510  DEN(l)  =  0. 

4520  DO  510  I-2rNP0INT 

4530  DE( I )  *  (E(I)-E(I-1))*DXI 

4540  DEP(I)  =  (EP< I )-EP(I~l ) >  *DXI 

4550  DEN(I)  *  (EH(I)-EM( 1-1 ) > *DXI 

4560  DENAX  *  MAX ( DE ( I ) r DENAX ) 

4570  DENIN  *  NIN(DE<  DrDEHIN) 

4580  DEPNAX  *  NAX  ( DEP  (Dr  DEPNAX ) 

4590  DENHIN  =  HIN(DEH< I) rDENNIN) 

4600  510  CONTINUE 

4610  PRINT r 'WANT  TABULATED  OUTPUT  ?* 

4620  READrANS 

4630  PRINT  692 r DENAX rDENINr DEPNAX rDENNIN 

4640  PRINTr 

4650  IF(ANS.NE.IHY)  GO  TO  530 

4660  PRINT  694 

4670  DO  520  I  “ISKIPrNPOINTrlSKIP 

4680  PRINT  630 r I r DE( I ) r DEP( I ) r DEH( I > 


4690  520  CONTINUE 

4700  PRINT  , 

4710  PRINT , 

4720  PRINT » ‘WANT  PLOTS  ?* 

4730  READ, ANS 

4740  IF(ANS.NE.IHY)  GO  TO  999 

4750  530  CONTINUE 

4760  PRINT, 'WANT  PLOT  OF  XX-DE  ?* 

4770  REAO,ANS 

4780  IF(ANS.EQ.IHY)  CALL  PL0T2!XX»DE»NP0INT> 

4790  PRINT, 'WANT  PLOT  OF  XX-DEP  ?* 

4800  READ, ANS 

4810  IF(ANS.EQ.IHY)  CALL  PL0T2(XX,DEP,NP0INT) 

4820  PRINT , 'WANT  PLOT  OF  XX-DEH  ?* 

4830  READ, ANS 

4840  IF(ANS.EQ.IHY)  CALL  PL0T2(XX,DEN,NP0INT) 

4850  PRINT , ‘WANT  REPLOT  ?* 

4860  READ, ANS 

4870  IF(ANS.EQ.IHY)  GO  TO  530 

4880  999  CONTINUE 

4890  CALL  PTIME(PTU) 

4900  PRINT  640, (PTU-PTI >{3600, 

4910  STOP  •  •  >  i. 

4920  600  FORMAT (4X» 1 1  * *5X» *EMAG*  »6X»  *EPHS* ,6X, ’ECOM'//) 

4930  610  FORMAT ( I6,3F10»&)  - 

4940  620  FORMAT  !5X»  *1 ' i'  .FMAGAM*,*  EPHSAV*,*  ECOMAV*//) 

4950  630  FORMAT (16, 3Flc£iK' 

4960  640  FORMAT ( *  PTU-SEC  =  *»Fl6.2> 

4970  650  FORMAT!//*BASE*:,5X,A60> 

4980  660  FORMAT! 'CASE  ■ ,I3,lX,A60//>  ' 

4990  665  FORMAT!/*  CASE',7110) 

5000  670  FORMAT! *  MAXC  *  *,7F10.3> 

5010  671  FORMAT! *  PEF  =  *,7F10.3> 

5020  672  FORMAT! *  EMMX  =  *,7F10,3) 

5030  673  FORMAT! *  EMMN  *  *,7F10.3> 

5040  674  FORMAT! *  EPMX  =  *,7F10.3) 

5050  675  FORMAT!*  EPMN  =  * ,7F10.3> 

5060  676  FORMAT! *  ECMX  =  *,7F10,3> 

5070  677  FORMAT!*  ECMN  =  *,7F10.3) 

5080  680  FORMAT ! /7X, *EMX* ,7X, *EMN* ,7X, 'EPMX* ,6X, * EMMN */2X,4Fl 0.3// 

5090  l  7X, *TE1 * ,7X, *TP1* ,7X, *TM1*/2X, 3F10,3) 

5100  682  FORMAT! *  TAR  =  *,7F10.3) 

5110  690  FORMAT !5X»  *  I  * ,6X, *E* ,9X, *EP* ,8X, 'EM* ) 

5120  692  FORMAT !/7X, *DEMX* ,6X, *DEMN* ,6X, ’DEPMX* ,5X, *DEMMN*/2X,4F10.3) 

5130  694  F0RMAT(5X, *  I* ,6X, *DE* ,8X, *DEP* ,7X, *DEM* ) 

5140  END 

5150  SUBROUTINE  READIN 

5160  PARAMETER  NC  *  10, N1  =  NC-1,NP  =  200 

5170  DIMENSION  TDUM120) ,C1 ! 4) ,C2!4) ,C3! 4) 

5180  CHARACTER  TITLE*60,TITL»20(3,NC> 

5190  CHARACTER  FILE*12,FMTF*9/9H!T12,1HJ )/,ANS*l 

5200  EQUIVALENCE  !TITLE,TITL) 

5210  COMMON  /INPUT/  XFINAL,DX,NTPLOT,ISKIP,NIBASE,NICOMP, 
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5220  S  ICNT  , NPOINT  , DT(NC) ,TAR(NC)  »NPTS(NC) . 

5230  C0HH0N  /ARRA1/  X(NP,NC) ,Y(NP,NC) »XX(NP> » YY(NP»NC) » 

5240  t  MAXC (N1 ) ,PEF  ( N1 ) , ENMX (N1 ) , EMMN (N1 ) »EPMX(N1 ) » 

5250  i  EPMN(N1 ) , ECMX (N1 ) , ECMN (N1 ) » TITLE (NC) 

5260  DATA  NOE, NAFT/0400000000000, 0403700000000/ 

5270  PRINT » ‘READ  XFINAL»DX»NTPLOT»SEARV»ISKIP»NIBASE»NICOHP* 

5280  READ,XFINAL,DX,NTPLOT,SEARV, ISKIP,NIBASE,NICOMP 

5290  ICNT  =  0 

5300  NPTM  =  NP 

5310  NPOINT  =  XFINAL/DX  +  1 

5320  IF(NPOINT.LE.NP)  GO  TO  5 

5330  XFINAL  =  (NP-1)  *  DX 

5340  PRINT  300, NP0INT»NP, XFINAL 

5350  NPOINT  =  NP 

5360  5  CONTINUE 

5370  PRINT » ‘NPOINT  =* f NPOINT 

5380  10  ICNT  =  ICNT  +  1 

5390  PRINT, 'READ  NSORCE.NFILE* 

5400  READ,NSORCE,NFILE 

5410  IF(NSORCE.EQ.O)  GO  TO  200 

5420  IF(NFILE.LT.l)  GO  TO  700 

5430  IF (NSORCE-1 )  200,20,100 

5440  20  REUIND  1 

5450  IF(NFILE.EQ.l)  GO  TO  40 

5460  DO  30  I  =1,2*(NFILE-1) 

5470  30  READ(1 »END=10) 

5430  40  READ(l)  NPTS< ICNT) ,DT( ICNT) , Cl ,C2,C3 

5490  NPT  =  MIN(NPTSdCNT), NPOINT) 

5500  IFtSEARV.LE.O.)  GO  TO  70 

5510  NPS  =  MIN(NPTS( ICNT) ,NP) 

5520  READ(l)  (XX(I) ,1=1, NPS) 

5530  DO  50  1=1, NPS 

5540  IF(XXd)-SEARV)  50,60,60 

5550  50  CONTINUE 

5560  60  CONTINUE 

5570  NSTRT  =  1-1 

5580  TARdCNT)  =  DT(ICNT)tNSTRT 

5590  NPT  =  NININPTS (ICNT) -NSTRT, NPT) 

5600  BACKSPACE  1 

5610  READd)  (SKIP,  K=1  ,NSTRT-1 ) ,  <Y(I,  ICNT),  1=1,  NPT) 

5620  GO  TO  80 

5630  70  CONTINUE 

5640  READd)  ( Yd , ICNT) ,  1  =  1 , NPT) 

5650  TARdCNT)  =  0. 

5660  80  CONTINUE 

5670  NPTSdCNT)  =  NPT 

5680  NPTM  =  MIN(NPT,NPTH> 

5690  CALL  BCDASC(C1  ,TITL(1 ,1  )!T)  ,20) 

5700  CALL  BCDASC(C2,TITL(2, ICNT) ,20) 

5710  CALL  BCDASC(C3,TITL(3, ICNT ) ,20) 

5720  GO  TO  10 

5730  100  REWIND  2 

5740  IF(NFILE.EQ.l)  GO  TO  140 


5750 

5760 

5770 

5780 

5790 

5800 

5810 

5820 

5830 

5840 

5850 

5860 

5870 

5880 

5890 

5900 

T710 

5920 

5930 

5940 

5950 

5960 

5970 

5980 

5990 

6000 

6010 

6020 

6030 

6040 

6050 

6060 

6070C 

6080 

6090 

6100 

6110 

6120 

6130 

6140 

6150 

6160 

6170 

6180 

6190 

6200 

6210 

6220 

6230 

6240 

6250 

6260 

6270 


HO  130  I  =  1  f2> (NFILE-1 ) 

130  READ(2fEND=10) 

140  READ(2)  NPTSdCNT)  fDT (ICNT) > TDUM 
NPT  =  NIN(  NPTSdCNT) fNPOINT ) 

IF (SEARV.LE.O. )  00  TO  170 
NFS  =  MIN(NPTS(ICNT)fNP) 

READ < 2 >  (XX( 1 > f 1=1 fNPS) 

DO  150  I=1fNPS 
IF(XXd)-SEARV)  150f  160f  160 
150  CONTINUE 
160  CONTINUE 
NSTRT  =  1-1 

TAR (ICNT)  =  DT(ICNT)$NSTRT 
NPT  =  N IN ( NPTS ( ICNT) -NSTRT f NPT) 

BACKSPACE  2 

READ (2)  <SKIPfK=1 fNSTRT-1 ) f ( Y( I f ICNT ) f 1=1 » NPT) 
GO  TO  180 
170  CONTINUE 

READI2)  ( Yd » ICNT )  f 1=1 » NPT > 

TAR (ICNT)  =  0. 

180  CONTINUE 

NPTSdCNT)  =  NPT 
NPTM  =  MIN(NPTfNPTM) 

CALL  BCDASC ( TDUMf TITLE (ICNT )f60> 

GO  TO  10 

200  ICNT  =  ICNT-1 
DO  500  K  =1 f ICNT 
DT(K)  =  PT(K)*1000, 

TAR(K)  *  TAR(K)*1000, 

DO  500  I  =1 » NPTS  <  K )  + 1 
X(I.K)  =  BT(K)  *  (1-1) 

500  CONTINUE 


IF (NPTM.GE.NF'OINT)  GO  TO  600 
NPOINT  =  NPTM 
XFINAL  =  (NPOINT-l )*DX 
PRINT  31  Of NPOINT » XFINAL 
600  CONTINUE 
RETURN 

700  CONTINUE 

PRINT f' INPUT  FILE  ?* 

READfFILE 

IF(FILE.EQ*1H  )  GO  TO  200 
CALL  DETACH(NSORCEff) 

ENCODE(FILEfFMTF) 

CALL  ATTACH(NSORCEfFILEfIfOfISTATf) 
IFdSTAT . EQ.NOE»OR.  ISTAT . EQ.NAFT)  GO  TO 
PRINTf’ISTAT  =  'fISTATf*  FILE  *  fFILE 
PRINT  96f ISTAT 
96  F0R(1AT(2Xf012) 

GO  TO  700 
98  CONTINUE 


PRINTf'READ  NFILE 


98 
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